Experimental  signature  of  programmable  quantum  annealing 
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Quantum  annealing  is  a  general  strategy  for  solving  difficult  optimization  problems  with  the  aid  of 
quantum  adiabatic  evolution  01-  Both  analytical  and  numerical  evidence  suggests  that  under  ide¬ 
alized,  closed  system  conditions,  quantum  annealing  can  outperform  classical  thermalization-based 
algorithms  such  as  simulated  annealing  HI-  Do  engineered  quantum  annealing  devices  effectively 
perform  classical  thermalization  when  coupled  to  a  decohering  thermal  environment?  To  address 
this  we  establish,  using  superconducting  flux  qubits  with  programmable  spin-spin  couplings,  an 
experimental  signature  which  is  consistent  with  quantum  annealing,  and  at  the  same  time  inconsis¬ 
tent  with  classical  thermalization,  in  spite  of  a  decoherence  timescale  which  is  orders  of  magnitude 
shorter  than  the  adiabatic  evolution  time.  This  suggests  that  programmable  quantum  devices, 
scalable  with  current  superconducting  technology,  implement  quantum  annealing  with  a  surprising 
robustness  against  noise  and  imperfections. 


Many  optimization  problems  can  be  naturally  ex¬ 
pressed  as  the  NP-hard  problem  of  finding  the  ground 
state,  or  minimum  energy  configuration,  of  an  Ising  spin 
glass  model  0, 0, 
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where  the  parameters  hj  and  Jjk  are,  respectively,  local 
fields  and  couplings.  The  operators  crj  are  Pauli  matrices 
which  assign  values  {±1}  to  spin  values  {"f , 4-} -  Two  al¬ 
gorithmic  approaches  designed  to  address  this  family  of 
problems  are  directly  inspired  by  different  physical  pro¬ 
cesses:  classical  simulated  annealing  (SA),  and  quantum 
annealing  (QA). 

SA  0  probabilistically  explores  the  spin  configuration 
space  by  taking  into  account  the  relative  configuration 
energies  and  a  time-dependent  (fictitious)  temperature. 
The  initial  temperature  is  high  relative  to  the  system  en¬ 
ergy  scale,  to  induce  thermal  fluctuations  which  prevent 
the  system  from  getting  trapped  in  local  minima.  As  the 
temperature  is  lowered,  the  simulation  is  driven  towards 
optimal  solutions,  represented  by  the  global  minima  of 
the  energy  function. 

In  QA  0,0  the  dynamics  are  driven  by  quantum , 
rather  than  thermal  fluctuations.  A  system  implement¬ 
ing  QA  is  described,  at  the  beginning  of  a  compu¬ 

tation,  by  a  transverse  magnetic  field 
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The  system  is  initialized,  at  low  temperature,  in  the 
ground  state  of  fftrans,  an  equal  superposition  of  all  2^ 
computational  basis  states,  the  quantum  analog  of  the 


initial  high-temperature  classical  state.  The  final  Hamil¬ 
tonian  of  the  computation  is  the  function  to  be  mini¬ 
mized,  -f/ising-  During  the  computation,  the  Hamiltonian 
is  evolved  smoothly  from  Htrans  to  Hising, 

H(t)  =  A(t)Htians  +  B(t)Hlsins,  t  e  [0,  T],  (3) 

where  the  “annealing  schedule”  satisfies  A(fS),  B(T)  >  0 
and  A(T)  =  B( 0)  =  0.  If  the  change  is  sufficiently  slow, 
the  adiabatic  theorem  of  quantum  mechanics  predicts 
that  the  system  will  remain  in  its  ground  state,  and  an 
optimal  solution  is  obtained  0!3-  Similar  transforma¬ 
tions  with  more  general  Hamiltonians  are  equivalent  in 
computational  power  (up  to  polynomial  overhead)  to  the 
standard  circuit  model  of  quantum  computation  |13l  Il4| , 
and  offer  at  least  a  quadratic  speed-up  over  any  classical 
SA  algorithm  fl5j. 

Realistically,  one  should  include  the  effects  of  cou¬ 
pling  to  a  thermal  environment,  i.e.,  consider  open  sys¬ 
tem  quantum  adiabatic  evolution  |16h21I|.  An  imple¬ 
mentation  of  open  system  QA  has  recently  been  re¬ 
ported  in  a  programmable  architecture  of  superconduct¬ 
ing  flux  qubits  [2  21-1251 .  and  applied  to  relatively  simple 
protein  folding  and  number  theory  problems  j2t|  |27I|. 
Although  quantum  tunneling  has  already  been  demon¬ 
strated  [25|,  the  decoherence  time  in  this  architecture 
can  be  three  orders  of  magnitude  faster  than  the  compu¬ 
tational  timescale,  due  in  part  to  the  constraints  imposed 
by  the  scalable  design.  In  the  circuit  model  of  quantum 
computation  this  relatively  short  decoherence  time  would 
imply,  without  quantum  error  correction  |28l[29j.  that  the 
system  dynamics  can  be  described  by  classical  laws  [Hj. 
In  the  context  of  open  system  QA,  this  might  lead  one  to 
believe  that  the  experimental  results  should  be  explained 
by  classical  thermalization,  or  that  in  essence  QA  has  ef¬ 
fectively  degraded  into  SA. 

Here  we  address  precisely  this  question:  are  the  dy- 
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FIG.  1:  Connectivity  graph  of  the  degenerate  Ising  Hamilto¬ 
nian  used  in  our  experiments.  The  four  spins  in  the  central 
square  are  the  “core  spins”  [the  first  four  in  Eq.  <[4j)] ,  the  four 
peripheral  spins  are  “ancillae  spins”  [the  last  four  in  Eq.  (|4|)]. 
As  depicted,  the  local  fields  hj  of  the  core  spins  have  value 
+  1,  the  local  fields  of  the  ancillae  spins  have  value  —1,  and 
all  couplings  Jjk  are  ferromagnetic  with  value  1. 

namics  in  open  system  QA  dominated  by  classical  ther- 
malization  with  respect  to  the  final  Hamiltonian,  as  in 
SA,  or  by  the  energy  spectrum  of  the  time-dependent 
quantum  Hamiltonian?  We  answer  this  by  studying 
an  eight-qubit  Hamiltonian  representing  a  simple  opti¬ 
mization  problem,  and  show  that  classical  thermalization 
and  QA  make  opposite  predictions  about  the  final  mea¬ 
surement  statistics.  Our  Ising  Hamiltonian,  depicted  in 
Fig.  [1]  has  a  17-fold  degenerate  ground  state 
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Sixteen  of  these  states  form  a  cluster  of  solutions  con¬ 
nected  by  single  spin- flips  of  the  ancillae  spins  [Eq.  (|4a[)] , 
while  the  seventeenth  ground  state  is  isolated  from  this 
cluster  in  the  sense  that  it  can  be  reached  only  after  at 
least  four  spin-flips  of  the  core  spins  [Eq.  (libl)].  As  we 
show  below,  classical  thermalization  predicts  that  the  iso¬ 
lated  solution  will  be  found  with  higher  probability  than 
any  of  the  cluster  solutions,  i.e.,  it  is  enhanced.  Fur¬ 
thermore,  after  an  initial  transient,  faster  thermalization 
corresponds  to  a  higher  probability  of  finding  the  isolated 
solution.  Open  system  QA  makes  the  exact  opposite  pre¬ 
diction:  after  an  initial  transient,  the  isolated  solution  is 
suppressed  relative  to  the  cluster,  and  faster  quantum 
dynamics  yields  higher  suppression  (lower  probability). 
Our  experimental  results  are  consistent  with  the  open 
system  QA  prediction  of  the  suppression  effect ,  and  in¬ 
consistent  with  classical  thermalization.  We  next  discuss 
these  opposite  effects,  starting  from  the  classical  case. 

Let  pi  denote  the  probability  of  state  i  in  the  clus¬ 
ter  (l4al) ,  and  ps  the  probability  of  the  isolated  state  (SB. 
The  thermalization  dynamics  are  dominated  by  single 
spin- flips  in  our  experiment  (see  Appendix).  The  prob¬ 


abilities  pi  are  all  close  because  states  in  the  cluster  are 
connected  by  single  spin  flips^  so  we  consider  the  average 
cluster  probability  pc  =  Enhancement  of 

the  isolated  state  means  that  ps  >  pc-  Our  SA  numer¬ 
ics  show  that  this  is  indeed  the  case  for  different  update 
rules  and  cooling  schedules,  throughout  the  thermaliza¬ 
tion  evolution  (see  Appendix).  To  explain  this,  note  that 
the  general  features  of  a  thermalization  process  are  de¬ 
termined  by  the  spectrum  of  iLising  and  by  the  combi¬ 
natorics  of  state  interconversion.  Each  of  the  17  degen¬ 
erate  ground  states  can  be  reached  from  any  other  state 
without  ever  raising  the  energy  via  a  sequence  of  single 
spin-flips,  so  that  SA  never  gets  trapped  in  local  minima 
(Appendix). 

The  SA  master  equation  and  the  classical  thermaliza¬ 
tion  prediction  ps  >  pc  can  be  derived  from  first  prin¬ 
ciples  from  an  adiabatic  quantum  master  equation  [2l|. 
Let  Hs(t)  and  Hsb  =  J2a  Aa  ®  Ba  denote  the  system 
and  system-bath  Hamiltonians.  The  Lindblad  equation 
is 

p  =  -i[Hs,p]  (5) 
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where  Labja  =  \a)(a\Aa\b)(b\,  wa6  =  Eb  -  Ea,  (|a)} 

the  instantaneous  eigenbasis  of  Hs  for  spin  vector  a,  and 
lap{w)  =  §°°00dTe'lu’T We  are  interested 
in  the  thermalization  process  in  which  the  density  oper¬ 
ator  is  diagonal  in  the  computational  basis  of  spin  vec¬ 
tors.  The  system-bath  coupling  Hamiltonian  then  has 

the  form  HSb  =  Ere{+,-,*}  Xql ®  ’  where 

(t 4"  =  (ox  ±  iav)/ 2.  We  denote  by  a+  (aj)  the  spin  vec¬ 
tor  resulting  from  flipping  the  jth  spin  up  (down).  From 
here  we  arrive  at  the  classical  master  equation  for  the 
populations  pa  =  /w 

N 

~  Eo-)Par  ~  fj{Ea  -  Ear)P(}j  , 

7=1 r=± 

(6) 

and  the  detailed  balance  condition  f(Ea  —  Ea±)  = 
exp [~P(E  ±  -  Ea)]f(E  ±  -  Ea).  Eq.  ©  is  the  mas- 
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ter  equation  that  we  used  in  our  SA  numerics.  It  can 
also  be  used  to  derive  the  classical  thermalization  pre¬ 
diction  ps  >  pc-  To  this  end,  it  can  be  seen  directly 
that  the  isolated  state  is  connected  via  single  spin-flips 
to  8  excited  states  with  energy  —4,  giving  the  rate  equa¬ 
tion  ps  =  8/(— 4)pe  —  8f(4)ps.  In  contrast,  all  states  in 
the  cluster  are  connected  via  single  spin-flips  to  at  most 
4  singly-excited  states;  the  other  four  spin-flips  connect 
between  other  states  in  the  cluster  and  hence  conserve 
the  energy.  Thus  pc  «  2  (/(-4)pe  -  f(4)pc).  Compar¬ 
ing,  we  conclude  that  population  feeds  into  the  isolated 
state  faster  than  into  the  cluster  and,  given  that  initially 
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FIG.  2:  The  time-dependent  gap  between  the  ground  state 
and  the  lowest  six  excited  states  in  the  relevant  region  of 
the  experimental  QA  evolution.  After  time  t  =  0.5 T  the 
highest  energy  level  shown  corresponds  to  the  isolated  state. 
The  inset  shows  the  transverse  field  magnitude  A(t)  and  Ising 
Hamiltonian  magnitude  B(t)  used  in  our  experiments,  during 
the  same  time  interval. 


Ps  >  Pc,  we  always  observe  ps  >  pc ■  Simultaneous  dou¬ 
ble  spin-flips  do  not  change  this  conclusion,  and  higher 
order  simultaneous  spin-flips  are  physically  less  likely.  A 
complete  derivation  is  given  in  the  Appendix. 

We  next  analyze  the  corresponding  predictions  of  QA. 
A  crucial  difference  with  respect  to  SA  is  that  now  the 
relevant  energy  spectrum  is  given  by  a  combination  of  the 
final  Ising  Hamiltonian  and  the  transverse  field.  Conse¬ 
quently,  as  shown  in  Fig.  [21  the  degeneracy  of  the  ground 
space  is  lifted  for  times  t  <T.  The  isolated  state  has  sup¬ 
port  only  on  the  highest  eigenstate  plotted  during  the  sec¬ 
ond  half  of  the  evolution.  Given  that  the  system  starts  in 
the  ground  state,  the  isolated  state  is  suppressed  by  the 
energy  gap,  until  this  gap  vanishes  at  the  end  of  the  evo¬ 
lution.  The  isolated  state  remains  suppressed  nonethe¬ 
less,  since  transitions  to  other  low  energy  states  require 
at  least  four  spins-flips.  The  transverse  field  term,  which 
drives  simultaneous  spin-flip  transitions,  is  small  at  large 
t.  If  the  four  spins-flips  are  not  simultaneous,  these  tran¬ 
sitions  involve  excited  states  with  much  higher  energy, 
and  are  suppressed.  This  predicted  QA  suppression  of 
the  isolated  state  is  confirmed  by  our  closed  and  open 
system  quantum  dynamical  simulations. 

Our  experiments  were  performed  using  the  D-Wave 
One  Rainier  chip  at  the  USC  Information  Sciences  In¬ 
stitute,  comprising  16  unit  cells  of  8  superconducting 
flux  qubits  each,  with  a  total  of  108  functional  qubits. 
The  couplings  are  programmable  superconducting  induc¬ 
tances.  The  qubits  and  unit  cell,  readout,  and  control 
have  been  described  in  detail  elsewhere  [221  [2rJ| .  The  ini¬ 
tial  energy  scale  for  the  transverse  field  is  10GHz  (the 
A  function  in  Fig  [21 ,  the  final  energy  scale  for  the  Ising 


FIG.  3:  Left:  schematic  depicting  the  maximal  connectivity 
graph  (Ah, 4)  of  the  qubits  inside  a  unit  cell.  Right:  an  em¬ 
bedding  of  -Hismg  from  Fig.  Q]  (right). 


Hamiltonian  (the  B  function)  is  5.3GHz,  about  fifteen 
times  the  experimental  temperature  of  17mK  ss  0.35GHz. 
To  gather  our  data,  we  ran  each  of  the  144  embeddings 
4000  times,  in  batches  of  1000  readouts,  resetting  all  the 
local  fields  and  couplers  after  each  batch. 

A  diagram  of  the  experimentally  achievable  coupling 
configurations  is  shown  in  Fig.  [3]  (left).  The  experimen¬ 
tal  results  are  shown  in  Fig.  01  The  key  finding  that 
is  immediately  apparent  is  that  the  isolated  state  is  ro¬ 
bustly  suppressed,  in  agreement  with  the  QA  but  not  the 
SA  prediction. 

Is  it  possible  that  suppression  has  an  explanation  other 
than  QA?  The  main  physical  argument  along  these  lines 
is  that  a  systematic  or  random  bias  due  to  experimental 
imperfections  breaks  the  17-fold  ground  state  degeneracy 
and  energetically  disfavors  the  isolated  state,  thus  lower¬ 
ing  ps  if  the  system  thermalizes.  We  proceed  to  examine 
this  and  the  robustness  of  the  suppression  effect. 

First,  note  that  spin  numbers  j  =  1,...,8  must  be 
assigned  to  the  flux  qubits  before  each  experimental 
run.  One  of  the  144  possible  such  “embeddings”  al¬ 
lowed  by  the  symmetries  of  the  Hamiltonian  and  the 
hardware  connectivity-graph  is  shown  in  Fig.  [3]  (right). 
Second,  note  that  spin-inversion  transformations  H(t)  — > 
<jjH(t)cTj  commute  with  F/trans,  and  simply  relabel  the 
spectrum  of  both  Rising  and  H(t ):  if  a  certain  spin  con¬ 
figuration  has  energy  A,  then  the  corresponding  spin 
configuration  with  the  jth  spin  flipped  has  the  same  en¬ 
ergy  E  under  crJiFisingi7j.  Spin-inversions  also  commute 
with  the  spin-flip  operations  of  classical  thermalization. 
Therefore  all  of  our  arguments  for  the  suppression  of  the 
isolated  state  in  QA  and  for  its  enhancement  in  classi¬ 
cal  thermalization  are  unchanged.  Using  spin  inversions 
we  can  check  that  the  suppression  effect  is  not  due  to  a 
perturbation  of  the  Hamiltonian  such  as  a  magnetic  field 
bias.  Indeed,  by  performing  a  spin-inversion  on  all  N 
spins  we  obtain  a  new  Ising  Hamiltonian  where  the  iso¬ 
lated  state  is  that  with  all  spins-up.  If  a  field  bias  sup¬ 
pressed  the  all  spins-down  state,  then  it  would  enhance 
the  all  spins-up  state.  Figure  0]  rules  this  out.  We  also 
tested  cases  with  only  antiferromagnetic  couplings,  and 
with  random  spin-inversions.  The  results  for  one  such 
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FIG.  4:  Statistical  box  plot  [3lJ  of  the  experimental  results  for  pa  (left  columns)  and  pc  (right  columns).  The  total  annealing 
time  was  T  =  hps  in  each  run.  In  each  column  the  bar  is  the  median,  the  box  corresponds  to  the  lower  and  upper  quartiles, 
respectively,  the  segment  contains  most  of  the  samples,  and  the  +’s  are  outliers.  Cells  1-4  are  physically  distinct  8-qubit  unit 
cells  on  the  chip.  No  statistically  significant  variation  is  seen  as  a  function  of  the  unit  cell  number.  “Inv”  is  the  case  where  each 
local  held  hj  is  flipped  to  —  hj.  In  this  case  the  isolated  state  corresponds  to  the  state  |TtTTtT1T)>  the  opposite  of  Eq.  (|4b)l. 
While  this  has  a  small  effect  for  cells  2  and  3,  in  all  cases  the  isolated  state  is  significantly  suppressed,  as  predicted  by  QA. 
This  establishes  that  suppression  of  the  isolated  state  is  not  due  to  a  global  magnetic  held  bias. 


random  inversion  example  are  shown  in  Fig.  [5]  In  all 
cases  we  found  agreement  with  the  QA  prediction,  but 
not  with  classical  thermalization. 

Robust  suppression  holds  even  at  the  level  of  individual 
embeddings  and  spin  inversions.  We  found  that  ps  < 
3%,  while  pc  >  6%  for  each  of  the  thousands  of  such 
cases  we  tested  (the  highest  median  for  the  experiment 
in  Fig  [4]  is  0.004).  Thus  suppression  survives  breaking 
of  the  ground  state  degeneracy,  which  certainly  occurs 
due  to  the  limited  precision  of  ~  5%  in  our  control  of 
{hj.  Jjk}-  The  suppression  effect  is  robust  because  it 
does  not  depend  on  the  exact  values  of  these  parameters, 
but  on  the  relatively  large  Hamming  distance  between 
the  isolated  state  and  the  cluster. 

Finally,  we  consider  the  effect  of  increasing  the  anneal¬ 
ing  time.  Open  quantum  and  classical  systems  converge 
towards  thermal  equilibrium.  Therefore  if  the  cause  of 
suppression  is  the  QA  spectrum,  longer  annealing  times 
will  result  in  ps  increasing ,  approaching  its  Gibbs  dis¬ 
tribution  value.  This  would  not  be  the  case  if  ps  were 
governed  by  the  spectrum  of  iTfemg.  In  Fig.  [S]  we  com¬ 
pare  a  numerical  simulation  of  open  system  QA,  using  an 
adiabatic  Markovian  master  equation  (2l[ ,  with  classical 
thermalization.  The  quantum  prediction  of  increasing  ps 
is  confirmed  experimentally,  as  shown  in  Fig.  [5] 

We  thus  arrive  at  our  main  conclusion:  signatures  of 
QA,  as  opposed  to  classical  thermalization,  persist  for 
timescales  that  are  much  longer  than  the  single-qubit 
decoherence  time  (from  5pis  to  20ms  vs  tens  of  ns)  in 
programmable  devices  available  with  present-day  super¬ 
conducting  technology.  Our  experimental  results  are  also 
consistent  with  numerical  methods  that  compute  quan¬ 


tum  statistics,  such  as  Path  Integral  Monte  Carlo  a.  Our 
study  focuses  on  demonstrating  a  non-classical  signa¬ 
ture  in  experimental  QA.  Different  methods  are  required 
to  address  the  question  of  experimental  computational 
speedups  of  open  system  QA  relative  to  optimal  classical 
algorithms. 
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Appendix  A:  First-principles  derivation  of  the 
master  equation 

Here  we  derive  the  master  equation  used  by  SA  from 
first  principles  within  the  open  quantum  systems  formal¬ 
ism.  This  motivates  classical  SA  as  a  model  for  a  system 
dominated  by  classical  thermalization  of  the  final  Ising 
Hamiltonian. 

Let  Hs{t)  be  the  time-dependent  system  Hamiltonian 
and  Hsb  =  '}2aAa®Ba  be  the  system-bath  Hamilto¬ 
nian.  We  have  previously  established  that  the  Lindblad 
equation  within  the  rotating  wave  approximation  has  the 


a  M.  Troyer  and  T.  Roennow,  private  communication. 
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FIG.  5:  Statistical  box  plot  of  the  probability  of  the  isolated  state  for  a  fixed  set  of  qubits,  with  different  annealing  times. 
In  this  plot  each  of  the  144  possible  embeddings  is  averaged  with  the  same  embedding  after  a  complete  spin  inversion.  This 
compensates  for  global  magnetic  field  biases  (which  can  be  seen  in  Fig.  [4]  cells  2  and  3).  The  Ising  Hamiltonian  for  this  data 
was  obtained  by  applying  a  random  spin  inversion  to  Hising  from  Fig.Q]  The  probability  of  the  isolated  state  increases  with  the 
annealing  time  T,  in  contrast  to  the  classical  thermalization  prediction.  While  in  classical  thermalization  an  initial  distribution 
concentrated  around  the  cluster  can  also  result  in  suppression  of  the  isolated  state,  this  suppression  is  highly  unlikely  to  persist 
after  a  random  inversion. 


form  HH 

p  =  -i[Hs,p]  (Al) 

+  lot*p(p>ab)  Labtpp]Jab’a  —  —  ^L^ab  aLabtp,  pj 


a/3  a^b 


+EE  7a*/3(0)  AQa>^/?Z/]^a  2  {^'aa,a^'bh>/3j  P} 

a/3  ab  L 


We  show  in  Sec.  [B]  that  for  a  bath  in  thermal  equilib¬ 
rium  at  inverse  temperature  (3 

7 a*/3(— w)  =  e_|3aJ7^Q*(w)  ,  (A4) 

where 

/OO 

dTelUT {Bp(T)Bl( 0))  .  (A5) 

-OO 


where 


Lab, a  =  |a)(a|Aa|6)(6| 

(A2a) 

Llb,a  =  \b)(b\AUa)(a\ 

(A2b) 

W ab  =  Ea  , 

(A2c) 

We  assume  that  the  system-bath  coupling  Hamiltonian 
has  the  form 

=  E  (A6) 

j=lre{±,2} 


{|a}}  is  the  instantaneous  eigenbasis  of  H$  (we  have 
suppressed  its  explicit  time-dependence)  for  spin  vector 
a  =  {ai, . . . ,  ajv},  where  as*  G  {t,  4-}  •  and 

/OO 

dreiulT  (B^  (r  )Bp  (0))  (A3) 

-OO 

is  the  Fourier  transform  of  the  bath  correlation  function. 
The  star  adornment  on  the  first  subscript  (a*)  is  a  re¬ 
minder  that  the  first  operator  in  the  bath  correlation 
function  is  Hermitian-transposed.  We  have  ignored  the 
Lamb  shift  in  Eq.  (IA1I)  since  for  a  time-dependent  Lind- 
blad  evolution  it  amounts  to  a  small  perturbation  of  the 
system  Hamiltonian.  We  used  this  form  of  the  master 
equation  for  our  quantum  open  system  numerical  simu¬ 
lations,  as  detailed  elsewhere  [2lj. 


where  a±  =  (ax  ±  iav)/2,  we  identify  |f)  with  |0)  and  |4_) 
with  1 1) ,  and  where  we  neglect  higher-order  interactions 
of  the  form  ctJ  ®  c (g)  B or  above.  Since  Hsb  is 
Hermitian  we  also  have  B^A  =  B^\  B^  =  BjZ\ 

g^*  =  gjT\  g^*  =  gjZ\  and  where  the  asterisk  denotes 
complex  conjugation.  In  the  computational  basis  of  spin 
vectors  {a},  we  introduce  the  notation 

\af)  =  af\a),  (A7) 

which  denotes  either  a  flipping  of  dj,  or  0  if  either  er d" 
acts  on  a,j  =t  or  aj  acts  on  aj  =4-.  Then 

(a\af\b)  =  (af)ab=Sa!b±  ,  (AS) 
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FIG.  6:  Probability  of  the  isolated  state  for  numerical  simula¬ 
tions  of  classical  thermalization  (Metropolis  update  rule)  and 
open  system  QA  as  a  function  of  the  total  annealing  time  T. 
“Ideal”  vs  “perturbed”  corresponds  to  simulations  for  i/iSmg 
without  and  with  a  perturbation  which  increases  the  energy  of 
the  isolated  state  (red).  In  classical  thermalization  p3  always 
decreases  with  T,  while  it  increases  for  QA  in  the  ideal  case. 
It  remains  almost  constant  for  QA  with  the  perturbed  Hamil¬ 
tonian.  Even  if  the  isolated  state  is  suppressed  energetically 
due  to  a  perturbation  of  Hiaing,  fast  classical  thermalization 
can  still  enhance  its  probability.  QA  with  the  ideal  Hamilto¬ 
nian  gives  the  best  qualitative  fit  to  the  experimental  data. 
System-bath  coupling  in  the  QA  simulation  corresponds  to  a 
decoherence  time  of  150ns. 


where  the  <5  function  is  defined  to  evaluate  to  zero  also 
when  <jj  annihilates  \b).  We  are  interested  in  classical 
thermalization,  in  which  the  density  operator  is  diagonal 
in  the  computational  basis  {a},  so  we  set  pab  =  0  for 
a  ^  b.  Equation  ED  then  gives  pab  =  0.  Using  indexes 
a  =  irij)  and  P  =  (s,  k)  in  Eq.  (1A1I).  where  r,s  £  {±,  zj 
and  j,k  £  [1 , . . . ,  iV] ,  and  taking  the  diagonal  (a\  ■  |a) 
matrix  element,  the  Lindblad  equation  becomes 


Paa=  9jr)*9{kS)x  (A9) 

(r,j),(s,k) 

T.  ”f(r,j)*(s,k)(wab)  ( crk)abPbb  ((crjV)ba 

b\b^a 

~  7(r,j)*(s,fe)  (wha)  ab  Paa  ( ak)ba  ' 


vanish,  and  using  Eq.  (IA8I).  we  are  left  with 

N 

P“  =  EE  \9jr)\2  (7(rj)*(rj)(wao-r)paj-r 
3= 1  r=± 

7(r,j)*(rb)  (^aJa)Pa^  j  (A10) 

where  we  denoted  pa  =  paa ,  the  probability  of  spin  con¬ 
figuration  a.  We  can  furthermore  identify 

P(d  a ))  =  (r,j)*(r,j)(uara)  (Alla) 

P{a~r  ->  a)  =  \gjr)\2J(rJ)*(r,j)^aaTr)  (AHb) 

as  the  transition  probabilities,  so  that  Eq.  (IA10I)  becomes 
the  rate  equation 


N 

Pa  =  EZ  P(a7r  aK"r  -  P(°  ^  aj)Pa  '  (Al2) 

j=i r=± 

This  can  be  further  simplified  using  the  KMS  con¬ 
dition.  Indeed,  note  that,  using  Ba(r)  =  ctj^t)  in 
Eqs.  (1A3I)  and  (1A5I).  we  have 


7(±j)*(±,j)M  =  7(T,i)(T,i)*H-  (A13) 

Using  this  along  with  to  ±  =  —  w  ±  [Eq.  (lA2ell  and 

Uj  aj 

Eq.  El,  we  have 

—  (3cj  ± 

7 (±,i)U±j)K±a)  =e  7(T,HUTd)Ka±)  ■  (A14) 

Therefore  Eq.  (1A11I)  yields 

P{a^af)=e  P  'aa±  \gf]  |27(Td)*(=Fd)  K„±)  (A15a) 
P(af  °)  =  l^(T)|27(Td)*frb)Ka±)  ■  (A15b) 


This,  together  with  gy  ’  =  gy  ,  gives  the  detailed  bal¬ 
ance  condition  for  thermalization  dynamics 


P(a  —>  af)  -p (E  ±-Ea) 
P{af  ~^a)=e 


f3{Ea  -  Ea±) 
fi(Eaf  -  Ea)  ’ 


(A16) 


where  we  introduced  transition  functions  fj(AE).  which 
we  identify  with  the  transition  probabilities  in  Eq.  (1A15I) . 

We  can  now  rewrite  Eq.  (IA12I)  as  the  classical  master 
equation  that  we  used  in  our  SA  numerical  simulations 


Pa  = 


N 


£  £  («% 


Ea)PaV 


fj  {Ea  Ea r  )pa^j  ■ 
(A17) 


Note  that  the  sum  in  Eq.  EB  involving  the  resonant  con¬ 
tribution  7c*/3(0)  vanishes,  since  the  terms  LaaippL'bb  a 
and  b  {E\a  aLbb,p,  p}  cancel  after  taking  the  diagonal 
matrix  element.  Moreover,  since  Eq.  m  involves  only 
off-diagonal  terms  ( b  ^  a),  contributions  due  to  uz  all 


Appendix  B:  Correlation  functions  and  the  KMS 
condition 

Here  we  derive  the  detailed  balance  condition  Eq.  (1A4I) 
from  first  principles.  Our  calculation  closely  follows 
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Ref.  Hi],  but  differs  in  that  it  applies  also  to  non- 
Hermitian  bath  operators. 

The  correlation  function  of  a  thermal  bath  is  as¬ 
sumed  to  satisfy  the  KMS  (Kubo-Martin-Schwinger)  con¬ 
dition  [22! 

(Bt(r)^(O))  =  <^(0)Bt(r  +  iP))  .  (Bl) 


This  expression  has  the  advantage  that  it  also  applies 
to  operators  which  are  not  trace  class.  For  trace  class 
operators  the  KMS  condition  can  be  derived  assuming 
that  the  bath  is  in  a  thermal  state,  ps  =  e~^HB ,  where 
Hb  is  the  bath  Hamiltonian.  In  this  case: 


(Bi(r)Bp(O))  =  Tr[pBU^(T,0)BiUB(T,0)Bp] 

=  -Tr  \Bpe~^-iT)HBBle-iTHB] 

z 

=  —  Tr[.B/ge^T+i|3)J?B  B^  e~^T+i^HB e~^HB] 

z 

=  T \\psBpUUr  +  ip,  0)4  Ub(t  +  iP,  0)] 

=  (^(0)Ht(r  +  iP))  ,  (B2) 


where  Ub  is  the  bath  unitary  evolution  operator.  Note 
that 


(Bi(r)Bp( 0))  =  {Bp(-t  -  *P)Bt  (0))  .  (B3) 


If  in  addition  the  correlation  function  is  analytic  in  the 
strip  between  r  =  —  i|3  and  r  =  0,  then  it  follows  that 
the  Fourier  transform  of  the  bath  correlation  function 
satisfies  the  detailed  balance  condition  Eq.  (1A4D  as  we 
show  next. 


the  closed  contour  encloses  no  poles  (the  correlation  func¬ 
tion  {Bp{r)B\{ 0)}  is  analytic  in  the  open  strip  (0  —  i(3) 
and  is  continuous  at  the  boundary  of  the  strip  |3J]),  so 
that 

£(...)  =  0  (B5) 

=  /(" 0+/(" -)  +  Ij- ’ 

where  (...)  is  the  integrand  of  Eq.  (IBdl).  and  the  integral 
f  is  the  same  as  in  Eq.  (IB4I) .  After  making  the  variable 
transformation  r  =  —x  —  <(3,  where  x  is  real,  we  have 

/  (...)  =  -eP“  r  e~iux {Bp (x) (0) )  .  (B6) 

J  <—  J  —  OO 


Assuming  that  (Ba (±00  —  z|3).Ba(0)}  =  0  (i.e.,  the  corre¬ 
lation  function  vanishes  at  infinite  time),  we  further  have 
(...)  =  (...)  =  0,  and  hence  we  find  the  result: 


dre^iBpi-T-imlm 


(B7) 


= 


/OO 

-OO 


which,  together  with  Eq.  (IB4I).  proves  Eq.  (IA4I). 


Appendix  C:  Spectrum  and  ground  states  of  the 
Ising  Hamiltonian 


FIG.  7:  Contour  used  in  our  proof  of  the  KMS  condition. 


We  compute  the  Fourier  transform: 


/OO 

dreluJT  (BI(t)B0(  0)) 

-00 

/OO 

-00 


(B4) 


The  spectrum  of  the  8-qubit  Ising  Hamiltonian  we  con¬ 
sider  in  the  main  text, 


can  be  analyzed  by  first  considering  the  spectrum  of  the 
Hamiltonian  coupling  a  single  ancilla  spin  to  a  core  spin, 
i.e., 


To  perform  this  integral  we  replace  it  with  a  contour  inte¬ 
gral  in  the  complex  plane,  j>c  drelulT {Bp(— r  —  i(3)Hj, (0)), 
with  the  contour  C  as  shown  in  Fig.  [7]  This  contour  inte¬ 
gral  vanishes  by  the  Cauchy-Goursat  theorem  [33j  since 


The  spectrum,  with  the  core  (ancilla)  spin  written  first 
(second)  is 


lit) 

in) 

lit) 

in) 


-i 

-i 

3 

-1 


(Cl) 


The  minimum  energy  is  —  1  whether  the  core  spin  is  up 
or  down.  It  is  important  to  note  that  if  the  core  spin  is 
up,  the  minimum  energy  is  —1  whether  the  ancilla  is  up 
or  down;  this  will  give  rise  to  a  16-fold  degeneracy  when 
we  account  for  all  spins  below. 

We  analyze  the  core  spins’  energies  by  first  taking  into 
account  only  their  couplings.  That  is,  we  analyze  the 
ferromagnetic  Hamiltonian 

1 

/  \ 

2  4 

\  / 

3 


Denoting  by  s  the  number  of  satisfied  couplings  (both 
spins  linked  by  the  coupling  have  the  same  sign) ,  the  en¬ 
ergy  is  4  —  2s,  where  s  £  {0,2,4}.  The  ground  states 
of  this  Hamiltonian  are  the  configurations  |  tttt)  and 
14444).  Since  Eq.  (IC1I)  shows  that  the  minimum  energy 
of  a  core-ancilla  pair  is  —1,  when  adding  the  low  energy 
configurations  of  the  couplings  to  the  ancillae  the  mini¬ 
mum  energy  is  —8.  It  also  follows  from  Eq.  (IC1I)  that  the 
ground  state  configurations  of  the  full  8-qubit  Hamilto¬ 
nian  are 


14,44-14411 )  (C2a) 

ittttmt) ,  (C2b) 

where  the  first  (last)  four  spins  are  the  core  (ancillae) 
spins,  and  ||)  means  that  the  spin  can  be  either  up  or 
down.  The  all-spins  down  case  (IC2al)  results  from  the  |44) 
configuration  in  Eq.  EH),  while  the  16-fold  degenerate 
case  (lC2bl)  results  from  the  degeneracy  of  |ft)  and|4t)- 
An  important  feature  of  the  energy  landscape  of  the 
8-qubit  Hamiltonian  is  that  it  does  not  have  any  local 
minima.  This  can  be  easily  proved  by  showing  that  a 
global  minimum  can  always  be  reached  from  any  state 
by  performing  a  sequence  of  single  spin  flips  and  never 
raising  the  energy.  To  see  this,  consider  an  arbitrary  state 
of  the  system.  We  can  first  flip  all  the  ancillae  spins  to 
|{)  which,  according  to  Eq.  EH,  can  be  done  without 
raising  the  energy  (independently  of  the  state  of  the  cor¬ 
responding  core  spin) .  Then  we  can  flip  the  core  spins  in 
order  to  satisfy  all  the  couplings  between  core  spins,  ei¬ 
ther  making  them  all  |4)  or  all|  {) ,  whichever  requires  the 
fewest  spin  flips.  Again,  according  to  Eq.  EH,  this  op¬ 
eration  will  not  raise  the  energy  of  the  core-ancilla  pair. 
Hence,  the  final  state  is  either  the  isolated  ground  state 
1 44414414) ,  or  the  state  |ttt14444)  that  belongs  to  the 
degenerate  cluster  of  ground  state  configurations. 


Appendix  D:  Simulated  annealing  and  classical 
thermalization 

Here  we  report  the  results  of  our  numerical  simula¬ 
tions  of  classical  thermalization  (or  SA),  using  the  master 
equation  (1A17I) .  In  the  plots  below  we  used  the  Metropo¬ 
lis  update  rule  for  the  transition  probability  P(a  — >  a±). 
Explicitly,  if  A E  =  Ea±  —  Ea  is  the  energy  difference  for 

the  update,  the  transition  probability  is 

u  l  .  min  (1,  exp(— /3A-E))  •  (Dl) 
#  of  spins 

We  have  also  tested  other  update  rules,  such  as 
Glauber’s  {35} ,  and  the  results  are  essentially  unchanged. 
The  result  are  also  essentially  unchanged  when  using  a 
transition  probability  pmin  (1,  exp(— (3AE))  for  positive 
p  (such  as  when  simulating  a  continuous  time  master 
equation).  We  do  find  a  dependence  on  the  choice  of 
annealing  schedule,  i.e. ,  the  functional  dependence  of  the 
temperature  on  the  number  of  steps.  Three  different  an¬ 
nealing  schedules  we  used  are  shown  in  Fig.  [H  and  the 
corresponding  SA  results  are  shown  in  Fig.  [9l  The  prob¬ 
ability  ps  of  the  isolated  state  is  always  above  the  average 
probability  pc  for  a  state  in  the  cluster. 

It  might  be  argued  that  thermalization  at  constant 
temperature  corresponds  most  closely  to  the  experimen¬ 
tal  situation,  given  that  the  experimental  system  remains 
at  an  almost  constant  17mK.  This  is  modeled  by  the  ex¬ 
ponential  annealing  schedule,  which  rapidly  converges  to 
a  nearly  constant  temperature,  as  can  be  seen  in  Fig.  [8] 
On  the  other  hand,  the  energy  scale  of  the  Ising  model 
changes  during  the  QA  evolution  (see  the  Fig.  2  insert 
in  the  main  text),  and  the  cooling  schedule  is  deter¬ 
mined  not  by  the  temperature  alone  but  rather  by  the 
ratio  between  the  energy  scale  and  the  temperature.  We 
also  show  ps  and  Pc  for  an  exponential  schedule  with 
Metropolis  updates  and  different  numbers  of  steps  in 

Fig.  [TUI 


Appendix  E:  Classical  master  equation  explanation 
for  the  enhancement  of  the  isolated  state 

We  now  explain  why,  as  seen  in  the  numerical  simu¬ 
lations  shown  in  Figs.  [9]  and  [TUI  the  probability  of  the 
isolated  state  never  exceeds  that  of  the  average  of  the  16 
cluster  ground  states,  i.e.,  why 

1  16 

Ps~Y6^Pi'  (E1) 

We  are  interested  in  sufficiently  slow  thermalization  pro¬ 
cesses  (relative  to  spin  flip  rates),  so  that  states  connected 
by  single  spin-flips  have  similar  populations.  The  cluster 
of  16  degenerate  ground  states  and  the  isolated  ground 
state  are  connected  via  a  plateau  of  excited  states  with 
energy  —4. 
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FIG.  8:  Temperature  as  a  function  of  annealing  step  n  for 
three  different  schedules:  exponential  T(n)  =  Tir™xp,  linear 
T(n)  =  Ti/(nrun  +  1),  and  logarithmic  T(n )  =  T;/(log(n  + 
1) nog  +  1),  where  T,  is  the  initial  temperature,  Tf  is  the  final 
temperature,  and  rexp  =  (T//T;)1/ritot,  rim  =  (T//T<  —  l)/ntot 
and  rlog  =  (Tf/Ti  —  l)/log(ntot  +  1),  where  ntot  is  the  total 
number  of  annealing  steps. 


FIG.  9:  Probabilities  from  SA  for  the  three  different  schedules 
shown  in  Fig.  [8j  The  probability  of  the  isolated  state  pa  is 
always  higher  than  the  average  cluster  state  probability  pc  ■ 


Let  us  first  derive  a  rate  equation  for  the  isolated  state. 
A  single  spin-flip  of  a  core  spin  in  the  isolated  state  raises 
its  energy  by  4,  since  it  violates  two  couplings  between 
the  core  spins  and  corresponds  to  a  transition  from  ||4-)  to 
|t|)  (where  the  second,  ancilla,  spin  is  unchanged),  which 
doesn’t  change  the  energy  according  to  Eq.  m-  Like¬ 
wise,  a  single  spin-flip  of  an  ancilla  spin  in  the  isolated 
state  violates  no  couplings  and  corresponds  to  a  tran¬ 
sition  from  ||t)  to  |tt)  (with  the  core  spin  unchanged), 
which  raises  the  energy  by  4  according  to  Eq.  (1C  1  ll .  There 
are  8  ways  this  can  happen  (4  core  and  4  ancilla  spins 
can  be  flipped) .  Since  this  accounts  for  all  the  single  spin 
transitions,  Eq.  (1A17D  yields  the  rate  equation 

Ps  =  8/(— 4)pe  -  8f(4)ps  ,  (E2) 

where  pe  is  the  population  of  the  excited  states  with 
energy  —4.  Here  we  are  assuming  that  the  spin  flip 


FIG.  10:  Probabilities  from  SA  for  varying  total  numbers  of 
steps.  We  used  the  Metropolis  update  rule  with  an  expo¬ 
nential  schedule.  The  lines  correspond  to  100  (dotted),  1000 
(dashed)  and  10000  (solid)  steps.  The  upper  curves  (blue) 
correspond  to  pa,  while  the  lower  curves  (red)  correspond  to 
pc-  The  inset  is  a  magnification  of  the  boxed  part.  The  sepa¬ 
ration  between  the  probabilities  of  the  isolated  state  and  the 
cluster  increases  as  the  temperature  decreases. 


rate  is  the  same  for  all  sites  [corresponding  to  assuming 
gP  =  in  Eq.  (1A6I)]. 

We  next  derive  the  rate  equation  for  the  cluster,  once 
again  accounting  only  for  single  spin  flips.  For  states  in 
the  cluster  the  core  spins  are  all  up,  and  ancilla-spin  flips 
are  energy-preserving  transitions  between  states  in  the 
cluster.  For  core-spin  flips  we  need  to  analyze  two  dif¬ 
ferent  situations.  The  first  is  a  configuration  in  a  ground 
state  where  the  core- ancilla  pair  starts  as  |tt)  and  the 
core  spin  flips,  so  the  state  becomes  |  ||).  This  vio¬ 
lates  two  couplings,  with  energy  cost  4,  and  according  to 
Eq.  (IC1I)  the  energy  difference  between  these  two  states 
is  4,  so  the  overall  result  is  an  excited  state  with  energy 
0.  The  second  is  a  configuration  in  a  ground  state  where 
the  core-ancilla  pair  starts  as  1 and  again  the  core 
spin  flips,  so  the  state  becomes  |||).  This  again  violates 
two  couplings,  with  energy  cost  4,  but  costs  no  energy 
according  to  Eq.  m,  so  the  overall  result  is  an  excited 
state  with  energy  —4. 

Thus,  a  state  with  l  ancillae  with  spin  |  and  4  —  l 
ancillae  with  spin  f  connects  (via  single  spin-flips)  to  l 
excited  states  with  energy  —4  and  4 —l  excited  states  with 
energy  0.  To  write  a  rate  equation  for  pc  =  X)i=i  Pi/ 16 
we  shall  assume  that  all  excited  states  with  energy  0  (—4) 
have  probability  p( 0)  (pe),  and  all  states  in  the  ground 
state  cluster  have  equal  probability  pc-  Summing  over 
the  number  l  of  ancilla  with  spin  |  for  each  cluster  state, 
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the  rate  equation  is 

16  4  /A\ 

^2  Pi  =  (  i  )  (Z/(“4)Pe  “  1/(4) PC 

2=1  Z  =  0  '  ' 

+  (4  -  Z)/(-8)  p(0)  -  (4  -  0/(8)  pc)  (E3a) 
=  32(/(-4)pe-/(4)pc 

+  /(-8)p(0)-/(8)pc),  (E3b) 

so  that 

Pc  =  2(/(-4)pe-/(4)pc+/(-8)p(0)-/(8)pc)  ■  (E4) 

For  most  temperatures  of  interest,  relative  to  the  en¬ 
ergy  scale  of  the  Ising  Hamiltonian,  the  dominant  transi¬ 
tions  are  those  between  the  cluster  and  states  with  energy 
-4.  Transitions  to  energy  0  states  are  suppressed  by  the 
high  energy  cost,  and  transitions  from  energy  0  states  to 
the  cluster  are  suppressed  by  the  low  occupancy  of  the  0 
energy  states.  Then 

PC  «2/(-4)pe-  2/(4)  pc  •  (E5) 

To  show  that  ps  >  pc ,  assume  that  this  is  indeed  the 
case  initially.  Then,  in  order  for  pc  to  become  larger 
than  pS7  they  must  first  become  equal  at  some  inverse 
annealing  temperature  (3b  ps(|3')  =  pcift')  =  Pg,  and 
it  suffices  to  check  that  this  implies  that  ps  grows  faster 
than  pc-  Subtracting  Eq.  (1E5I)  from  Eq.  (1E2D  yields 

Ps  —  PC  =  6  (/(— 4)  pe  -  /(4)  Pg) 

=  ,  (E6) 

\Pg  P(e^g)J 

where  in  the  second  line  we  used  Eq.  (IA16I).  Now,  be¬ 
cause  the  dynamical  SA  process  we  are  considering  pro¬ 
ceeds  via  cooling,  the  ratio  between  the  non-equilibrium 
excited  state  and  the  ground  state  probabilities  will  not 
be  lower  than  the  corresponding  thermal  equilibrium 
transition  ratio,  i.e.,  —  >  =  e-4^  .  Therefore,  as 

we  set  out  to  show, 

Ps-Pc>  0,  (E7) 


implying  that  at  all  times  ps  >  pc- 


Appendix  F:  Degenerate  perturbation  theory 
explanation  for  quantum  suppression  of  the  isolated 

state 

We  can  understand  the  splitting  of  the  degenerate 
ground  subspace  of  the  Ising  Hamiltonian  -f/iSmg  by  treat¬ 
ing  the  transverse  held  Htrans  =  —  Xq=i  aj  as  a  pertur¬ 
bation  of  the  Ising  Hamiltonian  Hising  (thus  treating  the 


QA  evolution  as  that  of  a  closed  system  evolving  back¬ 
ward  in  time).  According  to  standard  degenerate  per¬ 
turbation  theory,  the  perturbation  Pg  of  the  ground  sub¬ 
space  is  given  by  the  spectrum  of  the  projection  of  the 
perturbation  on  the  ground  subspace.  Denoting  by 

n0  =  ami)®8 + ]r  itttnmxtttnmi  (fi) 

the  projector  on  the  17-dimensional  ground  subspace,  we 
therefore  wish  to  understand  the  spectrum  of  the  opera¬ 
tor 


Pg  —  n0 


(F2) 


The  isolated  state  is  unconnected  via  single  spin  hips 
to  any  other  state  in  the  ground  subspace,  so  we  can  write 
this  operator  as  a  direct  sum  of  0  acting  on  the  isolated 
state  and  the  projection  on  the  space  Hq  =  H0  —  (I4-)  (4-l)<S>8 
of  the  cluster 


Pg  =  -0  ©  K 


(F3) 


While  <jx  acting  on  any  of  the  four  ancillae  connects  two 
cluster  ground  states,  ox  acting  on  any  core  spin  of  a  clus¬ 
ter  state  is  projected  away.  Therefore  the  perturbation 
is  given  by  the  operator 


Pg=~  0  © 


(F4) 


where  the  sum  is  over  the  four  ancillae  spins. 

Denoting  the  eigenbasis  of  ax  by  |±)  =  (|t)  ±  |/))/-\/2, 
with  respective  eigenvalues  ±1,  the  transverse  held  splits 
the  ground  space  of  i7ising  lowering  the  energy  of  |tttt 
+  +  ++),  and  the  four  permutations  of  |— )  in  the  ancillae 
spins  of  |tttt++H — )•  None  of  these  states  overlaps  with 
the  isolated  ground  state,  which  is  therefore  not  a  ground 
state  of  the  perturbed  Hamiltonian.  Furthermore,  after 
the  perturbation,  only  the  sixth  excited  state  overlaps 
with  the  isolated  state.  The  isolated  state  becomes  a 
ground  state  only  at  the  very  end  of  the  evolution  (with 
time  going  forward),  when  the  perturbation  has  vanished. 


Appendix  G:  The  quantum  Singular  Coupling  Limit 
does  not  agree  with  the  experimental  results 

Interestingly,  an  open  system  QA  master  equation  in 
the  singular  coupling  limit  (SCL)  yields  results  in  qual¬ 
itative  agreement  with  classical  thermalization,  and  op¬ 
posite  to  our  weak  coupling  limit  (WCL)  master  equation 
(EU).  Here,  following  Ref.  |32,  we  present  a  derivation  of 
the  SCL  master  equation. 

We  consider  a  Hamiltonian  of  the  form: 

H(t)  =  Hs(t)  +  +  e~2HB  ,  (Gl) 
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where  we  take  the  interaction  Hamiltonian  Hi  to  have 
the  form  A  ®  B,  where  the  system  (A)  and  bath  ( B ) 
operators  are  both  Hermitian.  The  formal  solution  in 
the  interaction  picture  generated  by  H$  and  Hb  is  given 
by: 


p(t)  =  p{ 0)  -  ie  1 


ds 


Hi(s),p(s) 


(G2) 


jo  L  J 

Plugging  this  solution  back  into  the  equation  of  motion 
and  taking  the  partial  trace  over  the  bath,  we  obtain: 


^Ps(t)  =  -£  2  J  dsTiB  (  Hi{t),  Hi{s),p(s)  )  , 


(G3) 

where  we  have  assumed  that  Tr [pbB]  =  ( B )  =  0.  Under 
the  standard  Markovian  assumption  that  p(t)  =  ps(t)  <8> 
Pb  and  under  a  change  of  coordinates  s  =  t  —  r,  we  can 
write: 


jtPs(t)  =  e  2  dr  [(A(t)p(t  —  r)A(t  —  r)  (G4) 
- A(t)A(t  -  t )p(t  -  t))  (B(t)B(O))  +  h.c.] 


Fig.  6  of  the  main  text.  Namely,  the  SCL  master  equa¬ 
tion  predicts  that  the  isolated  state  is  enhanced  relative 
to  the  cluster  of  states.  Since  the  bath  and  system-bath 
coupling  dominate  the  system  Hamiltonian  in  the  SCL, 
it  predicts  that  decoherence  takes  place  not  in  the  instan¬ 
taneous  energy  eigenbasis  of  the  system  but  in  the  com¬ 
putational  basis.  Furthermore,  the  resulting  Lamb  shift 
in  this  limit,  diagonal  in  the  computational  basis,  pref¬ 
erentially  lowers  the  energy  of  the  all-up  and  all-down 
states  relative  to  the  other  computational  states.  To¬ 
gether,  these  two  effects  cause  the  isolated  state  to  be 
more  populated  than  the  average  population  of  the  clus¬ 
ter  at  the  end  of  the  evolution,  in  contradiction  to  our 
experimental  findings.  The  agreement  we  find  between 
our  experimental  results  and  the  WCL  master  equation 
instead  supports  the  idea  that  decoherence  takes  place  in 
the  instantaneous  energy  eigenbasis  of  the  system  and/or 
that  the  Lamb  shift  does  not  dominate  the  system  Hamil¬ 
tonian. 


Appendix  H:  Entanglement 


where  A(t)  =  Us(t)AUg(t)  and  where  we  have  used  the 
homogeneity  of  the  bath  correlation  function  to  shift  its 
time-argument.  We  change  coordinates  r  =  e2r'  and  ob¬ 
serve  that  under  this  coordinate  change  (B(t)B(O))  is 
independent  of  e.  We  assume  that  this  bath  correlation 
function  decays  in  a  time  tb  that  is  sufficiently  fast,  such 
that  tb  -C  t/e2.  This  allows  us  to  approximate  the  in¬ 
tegral  by  sending  the  upper  limit  to  infinity.  We  also 
assume  that  tb  r'e2,  which  forces  the  correlation  time 
of  the  bath  to  zero,  hence  its  spectral  density  to  become 
flat,  and  hence — using  the  KMS  condition — amounts  to 
taking  the  infinite  temperature  limit.  Under  these  as¬ 
sumptions,  we  can  now  take  the  e  — >  0  limit,  yielding  the 
SCL  master  equation 

foPsit)  =  +  HLS,p(t)} 

+7(0 )^Ap(t)A-^{A2,p(t)}^j  ,(G5) 

where 

/OO 

dr' e~iu,T'  (B{t')B{Q))  ,  (G6) 

-OO 


where  L/ls  is  the  Lamb  shift  (renormalization  of  the  sys¬ 
tem  Hamiltonian)  and  where  V  denotes  the  Cauchy  prin¬ 
cipal  value.  Thus,  even  if  Hs  is  time-dependent,  we  re¬ 
cover  the  same  form  for  the  SCL  master  equation  as  in 
the  time-independent  case  [32j|.  This  SCL  master  equa¬ 
tion  yields  results  in  qualitative  agreement  with  classical 
thermalization,  and  opposite  to  the  WCL  presented  in 


Another  interesting  question  is  the  possible  presence  of 
entanglement  during  the  evolution.  Answering  this  ques¬ 
tion  experimentally  requires  measurements  to  be  per¬ 
formed  during  the  annealing,  a  capability  that  is  absent 
from  the  device  used  in  our  experiments.  However,  when 
we  compute  the  concurrence  of  the  states  obtained  from 
our  master  equation  (1A1I)  during  the  QA  evolution,  which 
is  consistent  with  the  statistics  of  the  measured  output, 
we  find  it  to  be  finite,  as  seen  in  Fig.  |TT]  This  suggests 
that  entanglement  is  being  generated  in  our  experiments. 


FIG.  11:  Concurrence  generated  during  our  QA  simulations, 
between  a  pair  of  ancillae  qubits  ( “external  edge” )  and  a  pair 
of  core  qubits  (“center  edge”).  We  show  the  concurrence  for 
the  ground  state  and  for  the  Gibbs  state  (“thermal”).  Our 
master  equation  Hi  for  the  time-dependent  density  matrix 
gives  concurrence  values  between  these  two  extremes,  that  de¬ 
pend  on  the  system-bath  coupling  strength  used  in  the  simu¬ 
lation. 


ffLS  =  -A- 


f 


du>^(u})V 
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